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REMARKS 

By the present amendment, the abstract have been amended in the manner 
suggested by the Examiner. Claim 6 has been amended to clarify features thereof, 
as discussed below. The specification does not require any correction since the 
proper spacing is provided between the words, as illustrated in US Publication 
2008/0231272. 

With respect to rejection of claims 1-8 under 35 USC 1 12, first paragraph, as 
failing to comply with the enablement requirement because the term "pseudo inverse 
matrix" is not described or defined in the specification, applicants submit that the 
definition of the term "pseudo inverse matrix" is well known, as disclosed, e.g., in 
Wikipedia (copy attached). That is to say, since the definition of this term is general 
and well known, there is no need to recite the definition of this term in the 
specification. Thus, claims 1-8 are considered to be in compliance with 35 USC 112, 
first paragraph. 

By the above amendment, claim 6 has been amended to clarify the features 
of the present invention, in particular to provide the structural relationship between 
the sliding window recited and rest of the elements or steps recited in claim 6. Claim 
6 has been amended to recite (in relevant part):"... the control unit... (3) employs a 
sliding window to obtain the plurality of images, and (4) performs the scanning of the 
k-space at intervals of n echoes and suppresses artifacts by a low-pass filter after 
image reconstruction". The limitations (3) and (4) are described by way of example 
on page 15, lines 9-16. For example, in the present invention, in order to suppress 
artifacts occurring on an edge of an image due to under-sampling, data to which the 
sliding window has not been applied and which is acquired by displacing scanning 
lines for fear scanning lines for adjoining frames may be superimposed on each 
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other is employed. The temporal filter (low-pass filter) is applied to the time- 
sequential data, whereby artifacts derived from under-sampling can be suppressed. 
Applicants submit that such features as now recited in claim 6 are in compliance with 
35 use 112, second paragraph, and patentably distinguish over the cited art as will 
become clear from the following discussion. 

With respect to the rejection of claims 6-8 under 35 USC 112, second 
paragraph, as noted above, independent claim 6 has been amended to be in 
compliance with 35 USC 112, second paragraph. Accordingly, dependent claims 7 
and 8 are considered in compliance with 35 USC 112, second paragraph. 

As to the rejection of claim 1-5 under 35 USC 102(b) as being anticipated by 
Tasaka et al. P.N. 6,515,477; the rejection of claims 6-8 under 35 USC 102(b) as 
being anticipated by Meyer et al. P.N. 5,485, 086, such rejections are traversed 
insofar as they are applicable to the present claims, and reconsideration and 
withdrawal of the rejections are respectfully requested. 

Independent claim 1 (taken as an example) recites (in relevant part):". ..the 
processing unit (1) divides the image echoes and reference echoes into a plurality of 
groups, (2) uses the reference echo and image echoes preceding and succeeding 
the reference echo to calculate an estimation coefficient , and (3) uses the estimation 
coefficient to estimate unmeasured echoes lying among the image echoes in the k- 
space" (emphasis added). These features are described, by way of example only, in 
lines 1-25, on page 11, in lines 1-2, on page 12, in lines 8-13, on page 13, in lines 5- 
11, on page 14 of the specification and illustrated in FIGS. 4-7. 

The present invention provides a magnetic resonance imaging technology for 
efficiently suppressing artifacts in radial scanning. In the present invention, the 
processing unit 104 illustrated in FIG. 3 measures seventy-two echoes. The 
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positions of echoes in the k-spaces are indicated with solid lines as illustrated in FIG. 
4. In FIG. 4, among echoes indicated with solid lines, an echo 213 indicated with a 
boldface line is a reference echo. The echo is used to calculate an estimation 
coefficient, whereby unmeasured echoes 212 indicated with dot lines are estimated. 
The sum total of echoes including estimated ones comes to one hundred and 
twenty-eight echoes. The procedure will be described in conjunction with the 
flowchart of FIG. 6. To begin with, one hundred and twenty-eight echoes including 
the unmeasured echoes 212 are divided into eight groups one of which includes 
seventeen echoes (step 401). As shown in FIG. 5, marginal echoes of each group 
are identical to marginal echoes of two adjoining groups. In FIG. 5, an echo 214 is 
included as a marginal echo in each of the first group 301 and second group 302. 
The reference echo 213 is measured to be included in the middle of each group. 
Thereafter, an estimation coefficient A = [a1 , a2] to be used to estimate the 
unmeasured echoes 212 included in each group is calculated according to the 
expression (1) below (step 402). 



where S equals [SI, S2] T (a T denotes a transposed matrix), S-denotes a pseudo 
inverse matrix, and SI and S2 (column vectors) are echoes adjacent to the 
reference echo R (row vector). Thereafter, the estimation coefficient A is used to 
estimate an unmeasured echo according to the expression (2) below (step 403). 



where Su denotes an unmeasured echo, and S' equals [S'1, S'2] T where S'1 and 
S'2 denote echoes adjacent to Su. By performing the above processing, each echo 
is equally divided into Np parts as indicated with dashed lines 303 in FIG. 5. FIG. 5 is 
concerned with a case where Np equals 3. When an echo is divided, a more highly 



A=RS- 



(1) 



Su=AS* 



(2) 
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precise estimate is obtained than the one obtained when one echo is used as it is. 
Finally, measured echoes and unmeasured echoes that are estimated as mentioned 
above are combined and gridded. Thereafter, inverse Fourier transform is performed 
in order to reconstruct an image (step 404). 

As mentioned above, according to the present invention, part of unmeasured 
echoes is measured as a reference. An estimation coefficient is calculated using 
echoes adjoining the reference, and the unmeasured echoes are estimated using the 
estimation coefficient. Consequently, since an echo measured as the reference is 
one of the unmeasured echoes, an imaging time hardly increases. Moreover, since 
the reference is used to calculate the estimation coefficient , unmeasured echoes can 
be highly and precisely estimated compared to when they are estimated by simple 
interpolation without using the reference. This results in a magnetic resonance 
imager capable of suppressing artifacts with little extension of the imaging time. 

On the other hand, Tasaka et al. discloses in FIGS. 9-17 a radial trajectory 

formation in sequence. Tasaka discloses in col. 6, lines 60-68: 

first, data collection is performed on a first trajectory 1 along the kx axis, as 
shown in FIG. 9. Next, data is collected on a second trajectory 2 along the ky 
axis, as shown in FIG. 10. The trajectories 1 and 2 are orthogonal to each 
other, and have an angular difference of tt/2. Thus, the k-space is divided into 
four quadrants. Then, data are collected respectively on trajectories 3 and 4 
each of which has a direction bisecting the angular difference tt/2 between the 
trajectories 1 and 2, as shown in FIG. 11. Thus, the angular difference 
between the trajectories becomes tt/4. 

Furthermore, Takasa discloses that when the data collection is performed, 
every trajectory has a different direction than the immediately preceding one by an 
angle of or near rr 12. Thus, time phases of motion of the imaged object 300 are 
approximately randomly distributed over the k-space. Therefore, regularity among 



the view data associated with the motion directivity is diluted . Thus, Tasaka 
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discloses a signal collection method by collecting magnetic resonance signals along 
a plurality of linear radial trajectories passing through center of k-space. However, 
Tasaka fails to disclose a processing unit that divides the image echoes and 
reference echoes into a plurality of groups. Tosaka also fails to disclose a 
processing unit that calculates an estimation coefficient by using the reference echo 
and image echo preceding and succeeding the reference echo. Tosaka also fails to 
disclose a processing unit that estimates unmeasured echoes lying among the 
image echoes in the k-space. 

Furthermore, Tasaka discloses making Cartesian (grid) data from radial data 
by interpolation method called gridding. The interpolation by gridding sums radial 
data with the weight response to distance from the interpolation (grid) point. Thus, 
the weight is predetermined and is not dependent on the image. In the present 
invention disclosed in claim 1, the estimation coefficient is calculated by using image 
echoes. Therefore, as it was stated above, the interpolation suitable for each image 
can be precisely performed as compared with interpolation method disclosed in 
Tasaka. Thus, Tasaka fails to disclose or teach:". ..the processing unit (1) divides 
the image echoes and reference echoes into a plurality of groups, (2) uses the 
reference echo and image echoes preceding and succeeding the reference echo to 
calculate an estimation coefficient , and (3) uses the estimation coefficient to estimate 
unmeasured echoes lying among the image echoes in the k-space" (emphasis 
added), as disclosed in independent claim 1 of the application. Thus, claims 1-5 are 
not anticipated by Tasaka, and should be considered allowable thereover. 

With respect to the rejection of claims 6-8, independent claim 6 (taken as an 
example) has been amended to recite (in relevant part):". ..the control unit (1) detects 
the nuclear magnetic resonance signal by radially scanning a k-space, (2) produces 
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a plurality of images, (3) employs a sliding window to obtain the plurality of images , 
and (4) performs the scanning of the k-space at intervals of n echoes and 
suppresses artifacts by a low-pass filter after image reconstruction " (emphasis 
added). These features are described, by way of example only, in lines 5-25, on 
page 16, in lines 1-15, on page 17 of the specification and illustrated in FIGS. 8-11. 
FIG. 10 shows the results of imaging performed according to the present invention. 
For example, FIG. 10(C) shows an image produced by scanning the k-space at 
intervals of three echoes, applying the temporal filter (low-pass filter), and 
suppressing artifacts derived from under-sampling. Not onlv streak artifacts 
attributable to a gap between adjoining data items but also artifacts appearing on an 
edge of an image are seen being suppressed . One example of the present invention 
has been described on the assumption that the k-space is scanned at intervals of 
three echoes. At interval of how many echoes the k-space should be scanned is 
determined using the expression (3) below. 



where N denotes the number of echoes N at intervals of which the k-space is 
scanned, and Ne denotes a magnification by which a frame rate is increased by 
sharing echoes. The expression (3) is drawn out as described below. Assuming that 
a frame rate is multiplied by Ne by sharing echoes, a gap between adjoining data 
items varies in cycles of Ne/F (where F denotes a frame rate). Accordingly, a streak 
artifact having a frequency F/Ne occurs. When the k-space is scanned at intervals of 
N echoes, the number of echoes to be scanned while a scanning line makes one 
turn is decreased to 1/(N+1). Consequently, the frequency of a change in the gap 
between adjoining data items is multiplied by (N+1). Accordingly, the frequency of a 
streak artifact comes to F/Ne x (N+1). The expression (3) is drawn out under the 



(N+1)=Ne/2 



(3) 
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condition of (F/Ne x (N+1) = F/2) that the frequency of a streak artifact should be 

high (Nyquist rate F/2), where Ne is magnification ratio of frame rate by sliding 

window. 

Meyer et al. discloses in FIGS. 6A and 6B the temporal filtering of the K-space 
signals for fluoroscopic imaging. In FIGS. 6A and 6B the moving window 40 moves 
with time to combine the latest three detected signals for the fluoroscopic image. 
Thus, in Meyer the contrast of the image can be improved by making the weight of 
the latest signal large . In contrast, in the invention recited in claim 6, a low path filter 
is used in the time direction of each frame image after image reconstruction . 

In the constitution which measures at intervals of N-echoes, the claimed 
invention, as it was stated above, can eliminate artifacts with the low-pass filter by 
shifting the frequency of the artifacts to high frequency region under the expression 
(N+1)=Ne/2. Thus, Meyer fails to disclose or teach: "...the control unit (1) detects 
the nuclear magnetic resonance signal by radially scanning a k-space, (2) produces 
a plurality of images, (3) employs a sliding window to obtain the plurality of images , 
and (4) performs the scanning of the k-space at intervals of n echoes and 
suppresses artifacts by a low-pass filter after image reconstruction " (emphasis 
added), as recited in independent claim 6 of the application. Thus, claims 6-8 are not 
anticipated by Meyer, and should be considered allowable thereover. 

In view of above amendments and remarks, applicants submit that claims 
present in this application should now be in condition for allowance and issuance of 
an action favorable nature is courteously solicited. 

If the Examiner believes that there are any other points which may be clarified 
or otherwise disposed of either by telephone discussion or by personal interview, the 
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Examiner is invited to contact Applicants' undersigned attorney at the number 

indicated below. 

To the extent necessary, Applicants petition for an extension of time under 37 
CFR 1.136. Please charge any shortage in fees due in connection with the filing of 
this paper, including extension of time fees, to the Antonelli, Terry, Stout & Kraus, 
LLP Deposit Account No. 01-2135 (Docket No. 520.46403X00), and please credit 
any excess fees to such deposit account. 

Respectfully submitted, 

ANTONELLI, TERRY, STOUT & KRAUS, LLP 



/Alan E Schiaveili/ 
Alan E Schiaveili 
Registration No. 32,087 



IS/AES 

1300 North Seventeenth Street. Suite 1800 
Arlington, Virginia 22209 
Telephone: (703) 312-6600 
Facsimile: (703) 312-6666 



2009 04/27 16:00 FAX 0335371624 Polaire IPC Tokyo USAN 8)004/017 

Generalized inverse 

From Wikipedia, the free encyclopedia 

(Redirected from Pseudoinverse) 
^'Pseudoijiverse" redirects here. For the Maore-Penrose pseudomverse. sometimes referred to as "the 
pseudoinverse", see Moore- Penrose pseudoinverse. 

In mathematics, a generaUzed inverse or pseudoinverse of a matrix A is a matrix that has some properties 
of the inverse matrix of A but not necessarily all of them. The term "the pseudoinverse" commonly means 
the Moore-Penrose pseudoinverse. 

The purpose of constructing a generalized inverse is to obtain a matrix that can serve as the inverse in some 
sense for a wider class of matrices than invertible ones. Typically » the genemlized inverse exists for an 
arbitrary matrix, and when a matrix has an inverse, then its inverse and the generalized inverse are the 
same, Some generalized inverses can be defined in any mathematical structure that Involves associative 
multiplication, that is, in a semigroup. 



Contents 

• 1 Types of generalized inveraes 

■ 2 See also 

■ 3 References 

■ 4 External links . 



Types of generalized inverses 

The various kinds of generalized inverses Include 

■ one-sided inverse, tha£ Is left inverse and right inverse 

■ Drazin inverse 

■ Group inverse 

*• Bott-Duffin inverse 

Moore-Penrose pseudoinverse 

See also 

■ Inverse element 

References 
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Subject Classification, MadiSciNet search (http://www.ams.org/m9ihsciaet/Bearch/publicati0ns html? 
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'Moore-Penrose pseodoinverse 

From Wikipedia, the free encyclopedia 

In mathematics, and in particular linear algebra, the pseudoin verse ^ of an m X 71 matrix A is a 

generalization of the inverse matrix. t-^^ More precisely, this article talks about the Moore-Penrose 

pseudoin verse, which was independently described by E. H. Moore^^^ in 1920 and Roger Eenroset^^ in 
1955. Earlier, Fredholm had introduced the concept of a pseudoinverse of integral operators in 1903. The 
term generalized inverse is sometimes used as a synonym for pseudoinverse. 

A common use of the pseudoinverse is to compute a 'best fit' (least squares) solution to a system of hnear 
equations that lacks a unique solution (see below under Applications) , The pseudoinverse is defined and 
unique for all matrices whose entries are real or comiplex numbers. It can be computed using the singular 
value decomposition. 
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Definition 

The pseudoinverse A of an m-by-n matrix A (whose entries can be real or complex numbers) is defined as 
the unique n-by-m matrix satisfying all of the following four criteria: •^^^'^'^^ 

1 . AA A = A (AA need not be the general identity matrix, but it maps all column vectors of A to 
themselves); 
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2, A'^ AA = A (A is a weak inverse for the multiplicative semigroup); 

3, (AA AA + (AA + is Hermitiau); and 

4, (A + A) * = A ^ A (A A is also Heraiitian). 



Here M is the hermitian transpose (also called conjugate transpose) of a matrix M. For matrices whose 
elements are real numbers instead of complex numbers, M = M^. 

Properties 

Proofs for some of these relations may be found on the proofs page. 

■ The pseudoinverse exists and is unique: for any matrix A, there is precisely one matrix A that 
satisfies the four properties of the definition. M If A has real entries, then so does A ^ ; if A has 
rational entries , then so does A . 

» If the matrix A is invertible, the pseudoinverse and the inverse coincide: A = A ^ ^ [5]:243 

■ The pseudoinverse of a zero inatrix is its transpose. 

■ The pseudoinverse of the pseudoinverse is the original matrix: (A ) = A.'^^l-^'^^ 

■ Pseudoinversion comriiutes with transposition, conjugation, and taking the conjugate transpose: ^^^-^^^ 

A = "A+j ■ 

« The pseudoinverse of a scalar multiple of A is the reciprocal multiple of A : 
(aA) + =:a'~U + forQ^ =^ 0, 



If A and B are such that the. product AB is- defined and 

■ A has orthonormal columns (A * A = /, unitarity is a special case), or 
« B has orthonormal rows (BB * = /), or 

■ A is of full column rank, and B is of full row rank, 



then 



(AB) + = i5 A + . 



The third case does not cover the first two; an orthonormal (including non- square) matrix must be of 
full rank, but otherwise there is no assumption made on the matrix it multiplies. 

AA is the orthogonal projector onto the range of A (the space spanned by the colunm vectors of A); 
A A is the orthogonal projector onto the range of A * , which agrees with the orthogonal 
complement of the kernel (null space) of A; (1 - A A) is the orthogonal projector onto the kernel 
(null space) of A S-'^^ 

The kernel of ^ + is the orthogonal complement of the range of A; the image of A is the orthogonal 
complement of the kernel of A. 

If the pseudoinverse of A * A is already known, it may be used to compute A : 
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A = (A * A) + A * . 

« Likewise, if {AA ' ) '*'.is already known: 

A + = A * ( A4 ) + . 
■ The pseudoinverse can be computed via a limiting process: 

^+ = Imi(/IM 4- Siy^A* = lira AUAA* + SI)'^ 

(see Tikhonov regularization) . These limits exist even if (AA * ) " ^ and (A * A) " ^ do not 
existJ'^^^^ea 

* > 

" In contrast to ordinaiy matrix inversion, the process of takuig pseudoinverses is not continuous: if the 
sequence (A^) converges to the matrix A (in the maximum norm or Frobenius norm, say), then (A„) 
need not converge to A - * 

Identities 

The following idehties can be used to cancel certain subexpressions or expand expressions involving 
pseudoinverses. Proofs for these properties can be found in the proofs subpage. 
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Special cases 

Orthonormal columus or rows 

If A has orthonormal columns (A * A = /) or orthonormal rows (AA * = /) , then A = A * . 

Full rank 

If the columns of A are linearly independent, then A * A is invertible. In this case, an explicit formula is; 

A + = (A * A) - U * . 
It follows that A is a left inverse of A: A A = /. 

r 

If the rows of A are linearly independent, then AA * is invertible. In this case, an explicit foimula is: 

-A. A. ' ^A-A. ' ^ t 

It follows that A is a right inverse of A: AA = /. 
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If both columns and rows are linearly independent (that is, for square regular matrices), the pseudoin verse 
is just the inverse: 



- 1 





Scalars and vectors 

It is also possible to define a pseudoinverse for scalars and vectors. This amounts to treating these as 
matrices. The pseudoinverse of a scalar x is zero if x is zero and the reciprocal of x otherwise: 

I 0, if.r = 0; 

otherwise. 

The pseudoinverse of the null vector is the transposed null vector. The pseudoinverse of a non-hull vector is 
the conjugate transposed vector divided by its squared magnitude: 

if T = 0; 
otlier.wise. 

For a proof, simply check that these definitions meet the defining criteria for the pseudoinverse. 

Circiilant matrices 

For a Circulant matrix C the singular value decomposition is given by the Fourier transform , that is the 
singular values are the Fourier coefficients. Let JPbe the Fourier matrix, then 

Block matrices 

Optimized approaches exist for calculating the pseudoinverse of block structm-ed matrices. 

Finding the pseudoinverse of ^ matrix 

Using regular inverses 

t 

Let k be the raxik of a X n matrix A. Then A can be decomposed as A = BC, where 5 is a 771, X Ar 
matrix and C is a X 71 niatrix. Then 

r 

^ + = C * (CC * )-'^iB* B)-^B* . 

If A has full row rank, so that k = m, then B can be chosen to be the identity matrix and the formula reduces 
to A = A (AA )~ \ Similarly, if A has full column rank (that is, /fc = n), we have A + = (A * A) ~ U * . 

Tlie QR method 

How^ever, computing the product AA ' or A A or its inverse explicitly is unnecessary and only causes 
additional rounding errors and computational cost. Consider first the case when A is full colunrn rank. Then 
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all that is needed is the Cholesky decomposition 

a''a = r*r, 

where R is an upper triangular matrix. Multiplication by the inverse is then done easily by solving a system 
with multiple right-hand-side, * 

by back substitution. To compute the Cholesky decomposition without forming A * A explicitly, use the QR 

decomposition .of A,A= QR where g is a unitary matrix, <2 * = * = /, and i? is upper triangular. 
Then 

A* A = {QBy{QB) = H'Q*QR == H*IB. = 7?.*/?,, 

so R is the Cholesky factor of A * A. The case of full row rank is obtained by swapping A and A * . 

The general case and the SVD method 

* ■ 

A computationally simpler and more accurate way to get the pseudoinverse is by using the singular value 
decoih.position.[^]t4][7] if ^ _ u^y * singular value decomposition of A, tlien A + = V2+C/*.Fora 
diagonal matrix such as 2,. we get the pseudoinverse by taking the reciprocal of each non-zero clement on 
the diagonal, and leaving the zeros in place. In numerical computation, only elements larger than some 
small tolerance are taken to be nonzero, and the others are replaced by zeros. For example, in the Matlab 
function pinv, the tolerance is taken to be r = 8-max(/n,n)«max(2), where e is the machine epsilon. 

The cost of this method is dominated by the cost of computing the SVD, which is several times higher than 
matrix-matrix multiplication, if a state-of-the 'art implementation (such as that of LAPACK) is used. 

The above procedure shows why taking the pseudoinverse is not a continuous operation: if the original 
matrix A has a singular value 0 (a diagonal entry of the matrix 2 above), then modifying A slightly may 
turn this zero into a tiny positive number, thereby affecting the pseudoinverse dramatically as we now have 
to take the reciprocal of a tiny number. 

Updating the pseudoinyerse 

For the cases where A has full row or column rank, and the inverse of the correlation matrix (AA * for A 
with full row rank or A A for full colimm rank) is already known, the pseudoinverse for matrices related 
to A can be computed by applying the Sherman-Morrison-Woodbury formula to update the inverse of the 
correlation matrix, which may need less work. In particular, if the related matrix differs from the original 

one by only a changed, ad.ded or deleted row or column, incremental algorithms'^] exist that exploit the 
relationship. 

Similarly, it is possible to update the Cholesky factor when a row or column is added, without creating the 
inverse of the correlation matrix explicitly. However, updating the pseudoinverse in the general rank- 
deficient case is much more complicated. t^i'i 9] 

R 

Software libraries 
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High quality implementations of SVD, QR, and back substitution are available in standard libraries, such as 
-LAPACK. Writing one's own implementation of SVD is a major programming project that requires a 
significant numerical expertise. In special circumstances, such as parallel computing or embedded 
computing, however, alternative implementations by QR or even the use of an explicit inverse might be 
preferable, and custom implementations may be unavoidable. 

Applications 

See also: Hat rnatrix 

The pseudoinverse provides a least squares solution to a system of linear equations. ^^^1 Given a system of 
Uhear equations 

Ax = 6, . 

in general we cannot always expect to find a vectors which will solve the system; even if there exists such 
a solution vector, then it may not be unique. We can however always ask for a vector x that brings Ax "as 
close as possible" to Z?, i.e. a vector that minimizes the Euclidean norm 

l\Ax - bf. 

If there are several such vectors jc. we could ask for the one among them with the smallest EucUdean norm. 
Thus formulated, the problem has a unique solution, given by the pseudoinverse: 



This description suggests the followmg geometric construction of the pseudoinverse of an mxn matrix A. 
To find A-^ b for given binW^, first project b orthogonally onto the range of A, finding a point pib) in the 
range. Thenform A-i({p(Z»)}), i.e. find those vectors in R" that A sends topib). This wiU be an affine 
subspace of R'' parallel to the kernel of A. The element of this subspace that has the smallest length (i.e. is 
closest to the origin) is the answer A'^ bwe are looking for. It can be found by taking an arbitrary member 

1 

of A" (•{>(£>)}) and projecting it orthogonally onto the orthogonal complement of the kernel of A. 
Using the pseudoinverse and a matrix norm, one can define a condition number for any inatrix: 

cond(/l)= 11.41111.4-^11 

A large condition number impUes that the problem of finding least-squares solutions to the corresponding 
system of linear equations is ill-conditioned in the sense that small errors in the entries of A can lead to 
huge errors in the entries of the solution.!^^^^ 

Generalization 

In order to solve more general least-squares problems, one could try to define Moore- Penrose 
pseudoinverses for all continuous linear operator A : ^2 between two Hilbert spaces Hi and H2, 

using the same four conditions as in our definition above. It turns out that not every continuous linear 

Operator has a continuous linear pseudo-inverse in this sense. t^^] Those that do are precisely the ones 
whose range is closed in 
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Selection of a Convolution Function for Fourier 

Inversion Using Gridding 

John h Jackson, Craig H. Meyer, Dwight G. Niahimura^ Member, IEEE, and Albert Macovski, Fellow, !BEE 



Absiract-^ln fields ranging from radio astronomy to magnetic 
resonance imaging, Fourier inversion of data not falimg on a 
Cariesiati grid has been a prbblem. As a result, multiple algo- 
rithm^ bave been created Tor reconstructing images Troni non- 
uniform frequency; samples. In the technique known as grid- 
ding» the data samples are weighted for sampling density and 
cdn'volved with « HnHe kernel, then resanipled on a grid pre- 
paratary to a fast Fourier transform. This paper compares ihe 
utility of several convohition functions, including one that out- 
performs the optimal" prolate spheroidal wove function In 
some situations. 



I. Introduction 

IN FIELDS ranging front radio astronomy to comput- 
erized tomography and magnetic resonance imaging, 
spatial frequency data is nsed to generate iniages. In order 
to take advantage of the great computational speed af- 
forded by the fast Fourier transform II), the data must lie 
on a Cartesian grid; However, because of hardware con- 
straints or practical considerations » this is not always fea- 
sible. Accordingly, many algorithms have been devel- 
oped for reconstruction from nonuniform samples. Some 
methods use various interpolation schemes, such as neer- 
est-neighbor, bilinear interpolation, and tmncated sine 
function, finite impulse response interpolators [2]-[4]. 
Other techniques include gradient descent methods [5], or 
reconstruction using coordinate transfomnation [6], For 
data on Bipolar grid, such as may be encountered in com- 
puterized tomography or diffraction tomography* Ultered 
back-projectioii [7] may be used for reconstruction. 

In this paper we consider the algorithm known as 
'*gridding.'* In its earliest form, as first used by radio 
astronomers, the spatial frequency plane was divided into 
a grid^ and the point in the center of each "cell" was then 
assigned a value equal t9 the sum of all of the data points 
falling within the grid [8] . Later improvements included 
the use of the average value of the data within the cell [9] , 
or weighting the sampling points based on the distance 
from array points, such as with a Gaussian function [10]. 
To account for variations in the spectral sampling density, 
a further modihcation of the Gaussian weighting method 
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normalizes the sum of the coefficients to be applied to the 
sampling points to unity, These cell summing, cell aver- 
aging, and Gaussian methods have been compared by 
Thompson and Bracewell [11], 

A.n overview of the gridding operation is given fay 
O^Sullivan [12], who shows that the optimal gridding 
method is convolution with a sine (sin irxf-Tf-x) function 
of infinite extent, followed by sampling- onto a Cartesian 
grid. Practical considerations require that the infinite sine 
function^ be replaced with a finite convolving function. 
This paper compaxes the artifact introduced into the image 
for various convolving functions of different sizes, in- 
cluding the Kaiser-Bessel window and the zero-order pro- 
late spheroidal wave function (PSWF). We also show a 
convolving fimction that improves upon the PSWF in 
some circumstances. 

II. Griddino Algorithm 

Consider a CWO-dimensionBl flinction m {x, y) with Fou- 
rier transform M{Uy v) given by 

'nOc, ^) exp [-2irz(i« + try)] dxdy, 

— 1» 

(1) 

and a sampling function S consisting of P two-dimen- 
sional delta functions at positions Uj, Vj^ 

p 

Siu, I/) S ^6(w — Uj,v - Vj). C2) 



The sampled Fourier data is given by 



(3) 



In gridding, the sampled data is convolved with a function 
C(u, V) [such as a Gaussian, a sine, or a small finite win- 
dow] and sampled onto a unit spaced grid, 

= {[Af • 5] ♦ C} • m (4) 

where * denotes • two-dimensional convolution, and the 
shah or comb function ni(M, v) is defined as a sum of 
equally spaced two-dimensional delta functions; 



i J 



(5) 
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The corresponding reconstructed image wises is given by 
the inverse Fourier transform of M^cs* 

mscsO^r y) = {[m{x, y) * sQc, y)] - cOc, y)} * mix, y). 

(6) 

Here we can observe the effect of 5(m, v) on tlie recon- 
structed image. First, the inverse Fourier transform of ihe 
original sampling function 5(u, v), which we refer to as 
s (jc. y) , aiFects the aliasing of mCf , y) at a level that cannot 
be recovered. Thus, as expected, if the function Af(w, v) 
is not sufUciently sampled, the aljaslng cannot be cor- 
rected via postprocessing. We can, however, make a cor- 
rection for a nohuniform sampling density in S(w, v) by 
introducing an area density function p{u, v). By defining 

pC", v) » S(u, V) * C{k, V), (7) 

areas that are oversamplcd will have a large area density , 
while areas that are undersampled will have a small area 
density. As discussed by Bracewel] and Thompson [13], 
the principal transfer function (PTF) is given by S{u, 
v)/ p (tt, t'), and theprincipal response pattern by the two- 
dimensional Fourier transform of the PTF. Introducing 
the area density function into reconstruction, we gen- 
erate the sampled, weighted, convolved, and sampled Af. 



Mi 



swcs 



.{[..[ 



S* C 



]] • ^) ■ 



m. 



(8) 



Note that the weighted and convolved Ms can be consid- 
ered a moving weighted average where the sum of the 
weighting coefficients as determined by C(«, v) has been 
norniaiized to one. The corresponding image is given by 

mswcsOcp y) = iim{x, y) * [s «"* {s • c)]} ■ c) * m 

(9) 

where + ^ refers to a deconvolution. 

We now consider the effect of C(tt, v) on the recon- 
sinicled image. As noted by O'Sullivan [12], the optimal 
convolution function is an infinite sine, but this function 
is computationally impractical. A finite convolving func- 
tion will contribute sidelobes, which will be aliased bade 
into the image by the shah ftniction. Also, any roUoif in 
the central lobe of cQc^ y) will show up as an attenuation 
towards the sides of the Image. This rollofif can be cor- 
rected by dividing by cix, y). This flattens the response 
across the central lobe, but. amplifies the effective ampli- 
tude of the aliasing sidelobes, as shown in Fig. 1 for a 
Hamming window. As is usually done, we will also limit 
the inxage to a finite region of interest containing only one 
** replication" of the object. This is represented mathe- 
matically by multiplying by a two-dimensional rect or 
boxcar function '^n{x, y) where 



p if w < a 

n(x, y) = i ^ ^ . 

v^O otherwise 



5 and |y| < 0.5 



(10) 




(a) 





Ft£. 1. (a) The Iog-«£iile inverse Pouiier transfom oC a 'Hamming win- 
dow. The doned lines repres^tii Ihe edges of the image, with energy in the 
inverse Fourier transform ouutde of rhe doited lines a1ia$in|; back into {he 
linage, (b) The convoluUon rolloST correctLon for the image, (c) Th€ effec- 
tive inveise Fourier tmosfonn after multiplying by the lolloff oonection> 



We call the generated image m*. 



m*{x, y) = mswcs 



'r-|(x, y) 
y) 



and we refer to the method of generating this image as the 
gridding algorithm. The process is illustrated in Fig» 2. 

III. C0NV01.UT10N Function Comparison 

The utility of any convolution function is deicnumed 
by che amplitude and placement of the aliasing sidelobes 
after the image has been corrected for rolloff near the 
edges. For example, gridding using a n(jc) function will 
suffer from large aliased sidelobes in the image, but none 
of the sidelobes will alias into the center of the image. If 
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S(u,v) 
Signal 



C(u,v) 




c(x,y) m(u,v) 



Image 




Fourier 
Transfonn 



We will consider the following convolving functions; 
I) two-term cosine, 

Of + (1 — OC) COS f ^ « 



Fig. 1, The gridding aigoriihm. The signal i» sampled by a ftii\ctton5(u. 
v). The resuUmg dala samples are convolved by thcfuncUon lO, which 
is also used lo generate ihe area sampling .densl^ weighting. The con- 
irolveddaea are sampled onto a Cartesian grid by muhiplying by m(.u, v), 
Finally, ihe invtrse Fourier transformed daia arc divided by c(JC, y) lo com- 
pensais for roUoif in the inverse Fourier transform ot C{u, v). 



the inverse Fourier transform of the convolation function 



2) thnee-ienn cosine, 
« + /J cos 

3) Gaussian, 



+ (1 - o: - ^) cos uy 



4) Kaiser-Bess&I, 



~^Zo[/3Vl - C2u/iVf]. and 



rolls off too much within the bounds of ihe image, the 5) prolate spheroidal wave function, see I15]-ri7J. 



aliasing into the sides of the image will be greatly ampli- 
fied when the roUoff is compensated. 

We will compaie onc-dimenslonal inverse Fourier 
tradsforms from separable functions, so C(u, v) — 
C{u)C{v), and c(;c, y} - c(jc) c(y). For convenience in 
working with scale factors, define the width of the 
desired image to be I unit. The inverse of this, 1 unit of 
spatial frequency, is Jhe spacing between sampling points 
in the shah function. 

The performance measure thai we will use is the rela- 
tive amount of aliasing energy, including the effect of the 
rolloff correction . The corresponding functional to be 
minirhized is ■ 



1 



(12) 



and thi^ is the measure that we will use for comparison. 

if all regions within the image are not of equal interest, 
the performance measure can include a spatially varying 
weight . function , such as Schwab's w(jc, y) [1 — 
(2j:)^]'*[1 - (2y)^]* (where a is a design paraijieter) [14], 
which is included in the funcdonal 



R = 



where'^ is the region of interest. Our equation (12) is 
similar to (13) where all regions within the image are of 
equal interest (so y^{x, y) 1), except that (12) includes 
the rolloff correction effect on the image, and we are con- 
sidering separable convolving functions, which simplijies 
the problem to one dimension, 



All functions are defined over \u\ giving each a 

width with csi, ^, and a as free design parameters. The 
well-known Hamming and Hanning windows are exam- 
plea of the two-term cosine function (with a — 0.54 and 
Of = 0.50, respectively), and the Blackman window is an 
.example of the three-term cosine function (with a =0.42 
and /3 == 0.50) [18]. For a given window width and de- 
sired bandwidth B the truncated zero-order.PSWF of the 
first kind contains the least amount of energy outside of 
the desired passband, i.e., it minimizes 



J\x\>B 
t) — oo 



dx 



(14) 



dx 



The PSWF is quite difficult to compute, bat the Kaiser- 
Bessel function [19] is a good approximation, and both it 
and its inverse transform are easily calculated. The func- 
tion itself is based on /q, the zeto-order modified Bessel 
function of the first kind, and the inverse transform is 
given by 

sin ^/tt^W^j^"'^^^^ 



For each of the functions, the parameters a, /3, and a 
were varied to detetmine the best possible performance, 
as measured by J in (12), at each function width. The 
resulting values of Jaresliown in Fig. 3 for the parameter 
values given in Table I. The relatively poor results are 
because of Ihe difficulty in generating a finite function 
whose inverse Fourier iransform goes in^mcdiat©ly from 
a near unity passband to a very low amplitude stopband, 
without allowing for a transition band. The consequence 
is relatively large errors from the aliasing sidelobes that 
appear near the edges of the image where the division by 
c{x)^ to make the passband uniform, also amplifies the 
corresponding portions of the sidelobes. 
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Fig! 3. A comparison of several convolving functions of v&T7ing widths, 
showing the relBtlvo amounl of ynwanted sidelobe energy ihat wiH alias 
into the Image Jf that convalvinf fUxici'ion is used in ihegrkSdmg algonthm. 
Note that when the. two-lcrm coaine function reaches a widlb of Toor, or 
when the three-term cosine function reaches a -width or six, the first side- 
lolie of the fuDciion's Invwjie Fourier irensfonn occurs within the image. 



TABLE I 

The Pakameter Vauuss tor Each Function Type that Provide thb 
1*5a5T Relative Aliased Bnerqy when GniDDiiqCi qnto a E^cular Grid 

Three-Term cos 



Window 


Two'Teizn cos 






Gaussian 


Kaiser-Bessel 


Width 


a 


a 


0 


cr 


8 


1.5 


0.7600 


0/8701 


0.231 1 


0.4241 


1.9980 


2.0 


0.7146 


0.a099 


0.31Q8 


0,4927 


2,3934 


2.J 


0.61S!» 


0.6932 


0.4176 


0.4839 


3.3BO0 


3.0 


0.5534 


0.5995 


0,4675 


0.5 063 


4.2054 


3.5 


0.5)85 


0.53S3 


0.4831 


0.5516 


4.9107 


4.0 




0.d99K 


0,4«9l 


0.5695 


5.7567 


4.5 




0-4653 


0,4972 


0.56S2 


6.6291 


5.0 




0.4463 


0.4985 


0.3974 


7.4302 



A simple solution to this problem, alluded to by 
0*Sullivaii [12], is to increase the image fieldrof-view by 
sampling onto a smaller grid, and then ignore the outer 
portion of ihe image. By doubling the image ficld-of-view, 
tJie central lobe of c (jc) can be three times as wide (since 
it can partially alias around both sides of the image and 
still not enter the region of Interest), which makes possi- 
ble much smaller amplitude sidelobes [16]. Additionally, 
cix) does not taper as much within the region of Interest, 
so less Folloff correction is needed. Conversely^ if the for- 
mer method inquired Fourier tranfiformatlon ofwu^NxN 
image, subsampling the data in the manner proposed will 
require transformation of a 2N X 2N image, in addition 
to computing Mswcs at four times as many points. The 
computation of Adg^Q^ at each point may be somewhat 
easier, however, since the subsampling method permits 
the use of a smaller couyolvlng hinction for any given 
error bound. 

In addition to the functions previously discussed, we 
have added a convolving function of our own design. The 
fimction was generaiecl in an iterative manner. Starting 
with any even function, the function was inverse Fourier 




Fig. 4. Authors' width 2.S convolving tHinction (a) and the carresponding 
inverse Fourier iraaafomi (b). Note the suppression of Ifae portions of Ihe 
.uansfonm that wUl alUa into the center rsglon of ihe imajgo. 




3 3.5 4 
Window Width 

Fig. S. A compadSQTi oT several convolving functions of varying widths, 
Ghowing the relative amount of unwanted sidclobe enerfiy tbac will alias 
Into the ceaier region of an image. Thia center region contains Ihe entire 
fcgioo of Interest if the sampled data was gridded onto a 2x subsampled 
grid. 



transformed and those portions of the transformed func- 
tion that would alias into the region of interest were set 
to zero. The function was then Fourier transformed back 
to the original domain and spatially bounded. This pro- 
cess wa5 then repeated many times. This is similar to the 
method used to generate the PSWF where an even func- 
tion is repeatedly low-pass filtered and spatially bounded. 
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TABLB 11 

Thr Pak.amkit.r Values i or Each Fl nction Tvpf. that Provide the 
L&Aur ReLATive auaski3 Enkko'i' when" GfiiODfva aura a 2 x 

Sua^iAMt'l-ltD Orio 



Three-Term cus 



Window 


Two-Term co.s 






Gan.'v.^ian 


KuLser-Be.'isei 


wmih 


or 






a 


0 


1.5 


0.5273 


0.4715 


0.49 17 


0.2120 


6.6875 


2.0 


0.5125 


0-4149 


0.4990 


0.2432 


9.1376 


2.5 


0.5076 


0.4OJI 


0.4996 


0.2691 




3.0 ■ 


0.5068 


0.3934 


0.4997 


0.2920 


I3,90ft6 


3.5 


0,505 1 


0.3897 


0.4999 


0.3145 


16.2734 


4.0 




0.3950 


0.5000 


0.5363 


18.SS47 


4.5 






0.5000 


0.3557 




S.O 




0.382.1 


0.-5000 


0.37.1? 







Fig. 6. Rcconsirucccd imagcn a numerical phwniom ciJn.»;i.Ming of four elHp.ses. .sbowing ihc cft'ecl Ihw ihe convoimioa fuoe- 
Uon hva on Ihe reiuliing Iniugc, Intoje {□) whs iTcnn.«mc!ed usinga wiOih ihree KafscT-Bessel window 10 = 4,20541. wUhout 
subsamplmg. Image (b) is [Tie .«mc as image U), hut scaled tmin 0 lo 5% lo show ihc crn.>r In Ihtf rewnsinieuon. Image <c> 
wii$ reconsimcied un'm^ a wldih three, t«0-ienn costne window (O.SOdS + 0.4932 cos lmt/1,5)] on a KUbjsainpJcd grid, i*nd Is 
scaled ihe same as image [b): linage (d) wiih rectinsimcwd using a width ihree. Kuiscr-BesseJ windrtw tji «= 13.9068) on a 



The desired PSWF is the eigenfunclion of this operation 
that has the largest eigenvalue, so ail undesired portions 
of the function will gradually decay away relalive to the 
desired function. In Fig. 4 our width 2.5 function is 
shown. For wider convolving functions, the sidelobes wiil 
be Increasingly lower, arid many more iterations wiU be 



required to genemie the function. The function shown is 
the resuU after I OOO 000 iteration.s. 

As seen in Fig. 5. the reduction in aliased energy is 
rather dramatic when the gridding is performed on a sub- 
sampled grid. The parameters a, p, and a for each of the 
function widths are given in Table II. 
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Fig. 6 shows the effect that the convolving function has 
on a reconstructed image. Reconstructions were done 
using A width three Kaiser-Bessel window on a regular 
grid, a width three two-terra cosine window on a subsam- 
pled grid, and a width three Kaiser-Bessel window on a 
subsaraplcd grid. Improvements in the reconstruction er- 
ror can be seen with each successive technique. 

IV, Conclusion 

The gridding algorithm can be impleniented with vir- 
tually no adverse effects from aliasing sidelobes of the 
convolving function. This requires that the convolved data 
is sampled finely enough to yield an image field-of-view 
that is larger than the actual region of interest. We have 
found that oyersampling by a factor of two in each direc- 
tion is sufficient to yield excellent sidelobe suppression. 

The selection of a convolving function requires two pri- 
mary considerations. First, the performance of the func- 
tion, and second, the computation time required to gen- 
erate the function. Although a (diacretized) version of a 
convolving function needs to.be computed only once, this 
can still be a limiting consideration if the generating com- 
putation time is extremely long, as is the case with the 
PSWF and the authors' function. Only slightly poorer re- 
suits, are achieved with the easily computable Kaiser-Bes* 
sel window. Use of the Kaiser-Bessel window requires 
selecting the free parameter, jS. The best choice of /S, in 
the sense of minimizing (12), is given in Table U for sev- 
eral window widths. 
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